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Abstract. A wavepacket model for a system of free pi- 
ons, which takes into account the full permutation sym- 
r — metry of the wavefunction and which is suitable for any 
^\ phase space parametrization is developed. The proper- 
ON ties of the resulting mixed ensembles and the two-particle 
^~^ correlation function are discussed. A physical interpre- 

(— I tation of the chaoticity A as localization of the pions in 

C^ the source is presented. 

^"5 Two techniques to generate test-particles, which sat- 
T-H isfy the probability densities of the wavepacket state, are 
C^ studied: 
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1. A Monte Carlo procedure in momentum space based 
on the standard Metropolis technique. 

2. A molecular dynamic procedure using Bohm's quan- 
tum theory of motion. 

In order to reduce the numerical complexity, the separa- 
tion of the wavefunction into momentum space clusters 
is discussed. In this context the influence of an unautho- 
rized factorization of the state, i. e. the omission of inter- 
ference terms, is investigated. It is shown that the corre- 
lation radius remains almost unefFected, but the chaotic- 
ity parameter decreases substantially. A similar effect 
is observed in systems with high multiplicities, where 
the omission of higher order corrections in the analy- 
sis of two-particle correlations causes a reduction of the 
chaoticity and the radius. 

The approximative treatment of the Coulomb in- 
teraction between pions and the source is investigated. 
The results suggest that Coulomb effects on the corre- 
lation radii are not symmetric for pion pairs of different 
charges. For (7r~,7r~) pairs the radius, integrated over 
the whole momentum spectrum, increases substantially, 
while for (7r"'",7r"'") pairs the radius remains almost un- 
changed. 



1 Introduction 

It is a well confirmed experimental fact [1, 2] that pi- 
ons, which are produced in heavy ion collisions, dis- 
play Bose-Einstein correlations. These correlations are a 
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consequence of the quantum mechanical interference in 
the corresponding symmetrical ra-particle wave function. 
They contain a large amount of information about the 
statistical properties of the momentum- and configura- 
tion space distribution of the system, and thus provide a 
method to probe the source geometry. In contrast to the 
situation found in astronomy, where Hanbury Brown and 
Twiss first measured the radii of stars by analyzing the 
correlation functions of photons [3], the pion sources in 
heavy ion collisions are of higher complexity. Computer 
simulations are necessary to interpret the pion correla- 
tion functions. In addition such simulations have to take 
into account the constraints imposed by the experiment. 
Adequate models therefore have to be versatile and flex- 
ible enough not only to allow for a proper treatment 
of the quantum effects, but also have to be based on a 
realistic description of the reaction and to include the 
experimental conditions. 

The strategy to investigate the correlation phenom- 
ena in heavy ion physics is displayed in the flow diagram 
Fig. 1: 

1) The goal of the source generator (Box 1) is to pro- 
duce the phase space distribution of the pions at freeze- 
out time, which is defined as the moment at which the 
strong interactions have ceased to exist. Since the fire- 
ball, which is the dominant source of pions, is not static 
but expands after being created in a state of highest 
density and temperature, and since the pions are proba- 
bly emitted at every instant during the expansion phase, 
the freeze-out time represents an individual property of 
each pion. Although there exist codes which use a para- 
metrization of the pion - phase space distribution (Box 
la), a more realistic description of the pion - emission 
process is realized in codes which are based on micro- 
scopical models like BUU or QMD [4, 5] and which do 
not provide parametrical expressions for the source (Box 
lb). 

2) After they are produced by the source generator, 
the pions still do not exhibit Bose-Einstein correlations, 
since the presently available codes do not include quan- 
tum statistics. The first step to obtain correlated pions is 
to construct a symmetrized wavefunction (Box 2). This 
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Fig. 1. The procedures needed to obtain effective correlation pa- 
rameter (Box 8) from source models (Box 1). 



can be done in different ways: The most common method 
is to use the two-particle plane wave state in momentum 
space (Box 2al) 



g-ipi-xi/ft g-ip2-X2/ft 
g-ipi-xj/ft g-ip2-xi/ft 



*(Pl,P2) 

(1) 

where pj denotes the momentum and Xj the coordinate 
of the i'th pion. As soon as the phase space density 
increases, higher order correlations begin to contribute 
to the statistics of the pions and a symmetrization of 
Tx'th order becomes the more adequate description of the 
quantum state (Box 2a2). 

3) Once the source properties are defined and the 
wavefunction is symmetrized, the 2-particle correlation 
function in momentum space 



Cl2(Pl,P2) 



'Pl2(Pl,P2) 
^(P1)^(P2) 



(2) 



can be constructed numerically. Here, 'P(p) is the single- 
particle momentum distribution, 'Pi2(pi, P2) is the prob- 
ability to find simultaneously two particles with mo- 
menta pi and p2. There exist a couple of numerical codes 
which are based on the plane wave approximation (Box 
2al) and additionally include approximative treatments 
of final state interactions via partial wave expansions [6]. 
4) Alternatively it is possible to produce correlated 
pion distributions by means of a quantum mechanical 
event generator. At present there exist two techniques 



for the eventwise generation of correlated pions: In the 
codes of the first type the pion pairs are weighted with 
the factor [7] 

1 + cos ((xj - Xj) • (pi - Pj)/h) , (3) 

which is the statistical weight for two particles described 
by plane waves as in Eq. (1). Additional weights occur 
in case of interactions. In Sect. 2.3 we will discuss to 
which degree these methods suffer from the inconsis- 
tencies inherent in the plane wave picture, which is a 
semiclassical approximation of the quantum mechanical 
state. Another problem is due to the fact that corre- 
lations of higher order than 2 are ignored. This is not 
true for the event generators of the second type, which 
where developed to produce events of high multiplicities 
and which employ quantum statistical methods [8] or ap- 
proximations for the 7x-particle state [9] (Box 2a2). These 
methods, however, are only applicable for the simplest 
source parametrizations (Box la) and can therefore not 
be combined with event generators based on microscop- 
ical models (Box lb). All presently available methods 
have in common that they exclude an explicit time de- 
velopment of the pionic system, since the description via 
plane waves or quantum statistics implies the existence 
of stationary states. 

5) The advantage of generated (hard) events as com- 
pared to a direct computation of the correlation function 
(Box 3) is that all analysis procedures, which originally 
were developed for the experimental data (dashed box), 
can also be applied to the simulated data. This includes 
detector simulations and tracking procedures and there- 
fore allows for the best comparability with the experi- 
ment. 

6) In general the comparison of the resulting corre- 
lation function Eq. (2) with the experimental function 
only tests the applicability of the model but does not 
yield values for the source parameters. 

7) Therefore, an analytical model is employed to ex- 
tract some characteristic source parameters. A simple 
parametrization of the pion source (Box 1) is the static 
fireball with a Gaussian density distribution, which we 
may call the standard pion source [9] : 



P(r) 



i^R') 



exp 



R^ 



(4) 



Note that R is not the rms radius, instead Rj \2 is the 
Gaussian cr in each direction, so that the rms radius of 
the source is given by ^JZ/2R. For such a source and 
when the particles are described as plane waves (Box 
2al) and chaotic emission is assumed, the corresponding 
correlation function can be evaluated analytically and 
yields [10] 



Ci2{q) 



exp 



2h? 



(5) 



where q = |p2 — Pi| is the magnitude of the momen- 
tum difference and the radius R is used as fit parameter. 
"Chaotic" means that the emission coordinates of the 
particles are independent, i. e. the two-particle correla- 
tion function in configuration space is a delta distrib- 
ution Ci2(ri — r2) = 6{ti — T2). This is a semiclassical 
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approximation in the sense that quantum interference ef- 
fects in configuration space are neglected. The picture of 
a chaotic source is a good approximation in astronomy, 
since the star radii are much larger than the coherence 
length of the emitted light, but it becomes inadequate 
in heavy ion physics, where the de Broglie wavelength 
of the pions is of the same order as the source dimen- 
sion. In order to compensate for this oversimplification 
and to obtain a better agreement with the experimental 
results, a second fit parameter A is introduced, which is 
called chaoticity parameter and which accounts for a pos- 
sible contribution of correlations in configuration space. 
The parameter A equals 1 in case of total chaoticity and 
equals in case of total coherence. We may call 



Ci2{q) 



A exp 



2h? 



(6) 



the standard correlation function. 

8) Finally, the effective radius iZeff and the chaoticity 
A are obtained by a fit of the analytical model Eq. (6) 
to the correlation function Eq. (2). 

The interpretation of the radius parameters extracted 
from two-pion correlations under conditions present in 
heavy-ion collisions at ultrarelativistic energies (Box 7 in 
the flow diagram) has been the subject of many recent 
publications, see for example [12, 13, 14, 17, 16]. 

This work is intended to close some of the gaps still 
present in today's numerical modeling of pion correla- 
tions (Box 2) and event generation (Box 4). In Sect. 2 
we present a quantum mechanical description of nonin- 
teracting many-pion systems in terms of wavepackets. A 
wavepacket model was already proposed by Sinyukov et 
al. [15], which, however, cannot be applied to pions with 
arbitrary momentum distribution, since it makes no dis- 
tinction between classical- and quantum mechanical en- 
sembles in momentum space. A more general treatment 
was presented by Padula et al. [16], who developed a gen- 
eralization of Pratt's formula [17] in order to describe 
wavepackets in Wigner space. These models, however, 
neglect the dispersion of the wavepackets, so that a cor- 
rect time dependent treatment in configuration space is 
not possible. 

We will present a most general description which is 
independent on any parametrization in configuration- 
and momentum space (Box lb) and includes ra-particle 
correlations (Box 2b, Sect. 2). This model introduces 
a new characteristic parameter of the system, the ini- 
tial wavepacket width (Tq. In Sect. 2.2.1 we discuss how 
this parameter effects the probability distributions of 
the pions in configuration- and momentum space. In 
Sect. 2.2.2 we discuss the correlation function of the 
wavepacket model and demonstrate how the finite width 
of the single-particle state accounts for partial coherence 
in a quite natural manner. A justification of the model 
in terms of physical arguments is offered in Sect. 2.3. 
In Sect. 3 a Monte Carlo procedure which allows us to 
produce events of correlated pions in full symmetry and 
independent of the source parametrization is presented. 
An approximative treatment of states with multiplicities 
higher than sa 10 is discussed in Sect. 3.2, in Sect. 3.3 



the influence of neglected interference terms on the 2- 
particle correlation function is investigated. In Sect. 4 
a molecular dynamic procedure which allows for an ex- 
plicit time development of the system in configuration- 
and momentum space is presented. First we demonstrate 
how Bohm's quantum theory of motion is suitable to ex- 
tend the rules for the time development of a quantum 
mechanical ensemble average into rules for single repre- 
sentatives of the ensemble. Then, in Sect. 4.2 an approx- 
imative treatment of Coulomb interaction between pions 
and their source is presented. 

Since our approach is of pure quantum mechanical 
nature, relativistic effects are not included, and a consis- 
tent way to extent the formalism so that it covers rela- 
tivistic quantum mechanics is at present not available. 



2 Generalized model 

We consider a set of n pions which may be produced by 
one of the usual source generators, e. g. BUU [4] or QMD 
[5], after freeze-out, so that no strong interactions occur 
at later times. The state of the system is defined by the 
phase space occupation 



{(Xi,Pi),(X2,P2),...,(X„,P„)} 



(7) 



Naturally, these pions do not display Bose-Einstein cor- 
relations, since no symmetrization was taken into ac- 
count in the course of the simulation. A large number 
of such states therefore build up a classical ensemble. 
The idea is first to symmetrize the system, construct 
a 7x-particle wavefunction W and then to generate a set 
of n test-particles whose phase space distribution coin- 
cides with the probability density of!?. The coordinates 
of these test-particles will be written in small letters, so 
that 



{(Xl, Pl), (X2, P2), . . . , (x„, p„)} 



(8) 



denotes a single representative of a quantum mechanical 
ensemble which is defined by the pure state W. 



2.1 Wavepacket representation 

Each of the n pions defines a Gaussian wavepacket with 
the following configuration space representation: 

^,(x, 0) = {2.al)-^l' exp (^ - ^^^^) ■ (9) 

It is centred in Xj and the corresponding probability 
density 

(Xi-x)2- 



|V-,|^ = (2,r^^)-^/^exp 



2al 



(10) 



has a width (Tq which is so far arbitrary. Its meaning will 
be examined in the following subsections. The solution 
of the Schrodinger equation yields: 

i^i[-yL,t) = [2i:s\t))-^l^x 

(P^-x(t)-g) (x(t)-X.-£^)- ^ ^^^^ 

h 4s(t)(7o 



exp 



(11) 
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Certainly, the wavepacket moves with the group velocity 
Vi/m, and, as a result of dispersion, has a time depen- 
dent width: 



s{t) 



ht 



2mai 



(12) 



The 7x-particle wavefunction is the symmetrized product 

!?(xi,X2,...,X„,t) = ^==y' V'<^(i)(xi,t)V!<^(2)(x2,t) 

Vra! — 

••• ■^<T(n){^X^,t) , (13) 

sometimes called permanent in contrast to its antisym- 
metrical counterpart, the determinant, which is required 
to describe systems of identical fermions. The sum in- 
cludes all n\ permutations {o"(l), o"(2), . . . , (T[n)}. On ac- 
count of simplicity we assume that all pions are emitted 
at the same time t = 0. Otherwise, we have to substitute 
s[t) — > s(tj) and Pjt — > Pjrj, where Tj is the time after 
production of the i'th pion. 

The momentum space representation of a single par- 
ticle is obtained by the Fourier-transformation 



(j)i {p,t) = h-^'^ I d^xi)i{yi,t) 



exp 



= A/i exp 
with 



h^ 



dpi 



Xj • dpi 



— zp • X 

h 
pH 

2m 




(14) 
(15) 

(16) 



and dpi = Pj — p. It displays the remarkable property 
of stationarity, i. e. there is no motion (since there is 
no final state interaction) or dispersion, only a rotating 
phase. The constant width is 



h 

2^ 



(17) 



The 7x-particle wavefunction has to be constructed in the 
same way as in configuration space: 



(l}a{n){Pn,t) ■ 



(18) 



In a first step, we have described the ra-particle distrib- 
ution (7) by a symmetrical ra-particle wavefunction W or 
#, respectively. From now on, the coordinates (Xj, Pj) do 
not play the role of particle coordinates. Instead, our task 
is to find a set of n test-particles, whose phase space dis- 
tribution (8) is defined by the probability densities |!?|^ 
(for the coordinates) and |#|^ (for the momenta), respec- 
tively. Before methods to obtain these particles are pre- 
sented, we first want to discuss the statistical properties 
of these test-particles, which can substantially deviate 
from the distributions of the classical ensemble (7). 



2.2 Ensemble properties of the test-particles 

We consider a simple physical situation: The configura- 
tion space distribution of the pion wavepackets is given 
by the standard source Eq. (4) and the momentum space 
distribution by a (classical) Maxwellian 

V 2mT / ^ ^ 



/(P) = (27rmT)-^/2 exp 



with a given temperature T. Of course, our description 
does not require to use a particular parametrization, but 
it will help us to develop an intuitive understanding of 
the model and provides us with the possibility to cross- 
check the results of the numerical codes. 



2.2.1 Single-particle distributions 

The momentum distribution of the test-particles is de- 
fined by the mixed state, which is obtained by integration 
of the single-particle (pure state) density |^(p)|^ over the 
(classical) density distribution Eq. (19): 



V{p) = (27rmT) 



-3/2 



/ \Hp) 



exp 



f^l." 



(27rmTeff)"^/^ exp 



(-. 



\2mT^ffJ 



\2mTj 



(20) 



The distribution remains to be of Maxwellian type, but 
the temperature has increased to 

„ „ ft' 



teff 



4:max 



T + T„ 



(21) 



The quantum temperature Tq is a consequence of the zero 
point energy the particles gain when they are localized 
in the configuration space by wavepackets. Alternatively, 
one may regard the increase of the temperature as the 
result of a delocalization in momentum space. 

The same calculation can be done for the configura- 
tion space distribution, yielding 

V{^) = {-kR^)-^'^ f |V-(x)|2 exp (^] d^X 

^2- 



R^ 



i^R'y 



-3/2 



exp 



R^ 



(22) 



with 



*(Pl,P2,...,Pn,<) = -^^ <f>a{l){Pl,t)<f>a{2){P2,t) R=./WT2^\ 



(23) 

Certainly, the test-particles occupy a larger source vol- 
ume than the wavepacket centers, a consequence of the 
finite wavepacket size. We conclude that the introduction 
of wavepackets leads to a blurring of the single-particle 
distributions: Localization of the test-particles in small 
wavepackets introduces a zero point energy which leads 
to higher temperatures, a consequence of Heisenberg's 
uncertainty relation. On the other hand, large wavepack- 
ets yield a better localization in momentum space and 
therefore do not disturb the momentum distribution, but 
then the information in configuration space is lost. Thus 
the optimal width (Tq has to be a compromise between 
these two extremes. 
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2.2.2 Two-particle correlation function 

The 2-particle correlation function is defined by Eq. (2), 
where the nominator is obtained by integrating the 2- 
particle pure state density over the (classical) phase 
space density distributions 



•'12 



J ^fl2 



p(Xi)p(X2)/(Pl)/(P2)|#(pi,P2)|' 

X d^Xi d^X2 d^Pi d^P2 . (24) 



2 Re{I}) (25) 



We have 

I^(P1,P2)|' = ^ (|<^1(P1)|'|<^2(P2) 
+ I<^2(P1)|'|<^1(P2)|^ 

with the interference term 

^ = <l>l{Pl) <l>2{P2)(l>*2{Pl) (l>l{P2) ■ (26) 

Notice that the interference term, in contrast to the 
single-particle probability densities, depends on momen- 
tum as well as on configuration space coordinates, which 
requires to integrate |#(pi, P2)|^ over the complete phase 
space. The norm A/12 is 



^fl2 



exp 



/ l*(Pi,P 
(Pi 



2)|'dVd'P2 



P2)' 



4,72 



exp 



-(Xi 



X2)2 



4<72 



A/12 is also a function of the centre coordinates and has 
to be included in the integration Eq. (24). Physically, the 
second term in A/12 denotes the degree of overlap in phase 
space. It aquires values between and 1, thus the integral 
Eq. (24) can be solved in terms of a Taylor expansion of 
the inverse norm A/j^ and integrated term by term. The 
resulting expressions are quite cumbersome and provide 
little physical insight. Therefore, in subsequent sections 
we will treat the integration numerically, but beforehand 
we may get some rough impression about the solution by 
means of approximations. 

First we may examine the O'th order Taylor contri- 
bution, i. e. the result of Eq. (24) in the approximation 
A/12 = 1- We obtain 

Ci2(pi, P2) = 1 + exp f -^e«-(P^-PO^ ) (28) 



with the effective radius 



R. 



eff ■ 



(29) 



R^ + 2al—. 

Jeff 

Two points are remarkable: 

1. Like the standard correlation function Eq. (6), C12 
only depends on the absolute momentum difference 
|P2 — Pi I, but the chaoticity parameter A is fixed to 
1. As we will see below, this changes if Taylor terms 
of higher order are taken into account. 

2. The effective radius obtained by interferometry is not 
equal to the radius R of the single-particle distribu- 
tion Eq. (22). Instead, there is an additional contri- 
bution from the quantum temperature. We will see 
in Sect. 3 that higher order Taylor terms increase the 
value of -Reff- 



Finally we demonstrate how the finite size of the 
wavepackets effects the chaoticity A for our model. A 
more general treatment of partial coherence is given in 
[15]. For our argumentation we adopt a quantum statis- 
tical point of view and assume that due to the overlap 
of wavepackets there exist a finite number of K config- 
uration space cells. Different cells represent independent 
states. Pions which are emitted from identical cells rep- 
resent the coherent part of the correlation function and 
do not contribute to the chaoticity A, while pions from 
different cells are independent and increase the chaotic- 
ity. Therefore we define 



Number of coherent pairs 
Number of all pairs 



(30) 



Employing the fact that the statistical weight of a cell, 
which is already occupied by N bosons, is equal to N + 1, 
A can be evaluated by means of combinatorial methods 
as shown in the appendix. We obtain 



(31) 



K + 1 ' 



independent of the multiplicity of the system. In order 
to estimate K, we note that A, following Eq. (6), can be 
expressed as 



(27) A = Ci2(g = 0) - 1 . 



(32) 



We use the two-particle probability density Eq. (24) and 
assume that not only the momenta of the test-particles 
are equal, i.e. pi = p2, but that in addition also the 
wavepacket momenta are equal, i. e. Pi = P2. As a con- 
sequence the norm A/12 Eq. (27) becomes independent 
of the momentum and the integration can be carried out 
easily. We obtain 



2E 

j=i 



(-ly 



(^ + ^) 



3/2 



(33) 



It is not surprising that A, and hence the number of 
configuration space cells K, is a function of the ratio 
2^, i. e. the source dimension devided by the single-state 
dimension. For smaller wavepackets, the number of cells 
increases, the source becomes more incoherent and the 
chaoticity parameter A approaches 1. 

The approximation Pi = P2 is responsible for the 
fact that A does not depend on the momentum coordi- 
nates. A more accurate treatment would yield that A is 
not independent on the momentum space distribution, 
even in case of static sources. We will demonstrate in 
Sect. 3 that Eq. (33) nevertheless provides a reasonable 
estimate for the chaoticity. 



2.3 Physical legitimation of the generalized model 

The wavepacket model seems to introduce a new para- 
meter (7o into the physics of pion interferometry. Such a 
step has to be justified by physical arguments, i. e. it has 
to be demonstrated that the expansion of the formal- 
ism not only increases the complexity of the formula but 
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also allows to describe additional observable phenomena 
which otherwise remain hidden. 

In order to demonstrate the conceptual difference be- 
tween the usual picture which leads to the standard cor- 
relation function Eq. (6) and the wavepacket model, we 
write, using the Fourier expansion Eq. (14), the two- 
particle state Eq. (25) in the form 

l^(pi,P2)|'=fe-'x 

I / fd^xid^X2 (V'i(xi)e-'P-"'/'^V>2(x2)e-'P-"^/'^ 



+ V-2(xi) e-'P-"'/'^ V-i(x2) e-'P-"^/'^) 



(34) 



In the next step, we replace the configuration space 
states by delta distributions: 



), 



V'i(xj) — > 6{Xi - X, 
which yields 

l^(pi,P2)|'=fe-' 



g-ipi-Xi/ft g-ips-Xj/ft 



g-ipi-Xj/ftg-ips-Xi/ft 



(35) 



.(36) 



Now we have reached the plane wave picture in momen- 
tum space Eq. (1). After substituting this expression into 
Eq. (24) it is immediately obvious that the resulting ex- 
pression is identical with the (incoherent) nominator of 
the correlation function [10], leading to Eq. (5). The sub- 
stitution (35) therefore orthogonalizes the states in con- 
figuration space and also decouples them from momen- 
tum space, since ■0j(xj) contains momenta which are lost 
in 5(Xj — Xj). This procedure introduces a chaotic source 
(chaoticity A = 1). On the other hand, in the same ap- 
proximation the single-particle momentum distribution 
Eq. (20) would become the unity distribution. This is, 
as well as the obtained total chaoticity, in disagreement 
with the experimental results. To overcome these diffi- 
culties, the chaoticity A is introduced as a free parameter 
to account for partial coherent sources as a consequence 
of nonorthogonal states in configuration space. This at 
least in principle mimics the necessity to localize the par- 
ticles in momentum space. In the usual treatment of pion 
correlations, however, the picture of localized particles is 
adopted and simultaneously the pions are forced to sat- 
isfy a certain energy distribution, which is a violation of 
Heisenberg's uncertainty principle and must be regarded 
as a semiclassical approximation of the physical reality. 

The wavepacket picture, on the other hand, accounts 
for all quantum mechanical properties of the system in 
a quite transparent manner: The localization of the par- 
ticles in configuration space leads to a delocalization in 
momentum space and vice versa. The chaoticity A is not 
a free parameter, but a function of other system parame- 
ters, especially of the geometrical dimensions in config- 
uration space Eq. (33). From this point of view, (Tq does 
not represent an additional new parameter, but instead 
a parametrization of A. 

It is an interesting question whether any physical 
meaning can be attributed to (Tq. In the opinion of the 
authors the concept of wavepackets is a natural conse- 
quence of the processes that occur in the fireball: A pion. 



which is created in a certain area of the fireball, has a 
small probability to migrate over large distances in the 
fireball. More specifically: Since there exists a mean free 
path I for the pion to travel inside the fireball, it is ef- 
fectively localized within a certain volume at freeze-out. 
The mean free path I depends on the momentum, for 
200 MeV/c pions and in nuclear matter with density po 
it is expected to be Z = 1.3 fm [18]. On the other hand, 
the density of the fireball at maximum compression is 
around 3po and near freeze-out it is found to be 1/3 po 
[19]. Hence the mean free path I for 200 MeV/c pions 
in the fireball can be expected to be anywhere between 
0.4 fm and 4 fm. Assuming that two pions interfere over 
a distance sa I and move independently in case of much 
larger distances, I may roughly set the size of the local- 
ization volume and therefore of the wavepackets. There- 
fore, the measurement of the chaoticity A in connection 
with the source radius may, at least in principle, provide 
some information about the mean free path and hence 
the cross section for pion-absorption in nuclear matter. 



3 Monte Carlo method 

The procedure that will be presented in Sect. 3.1 is a 
generalization of Zajc's method [9] applied to the model 
of Sect. 2. The ra-particle state is described in momen- 
tum space. The high complexity of the permanent Eq. 
(18) requires methods which are different from Zajc's 
approximation, since the latter one is only applicable in 
case of real valued wave functions. Instead we exploit 
the fact that as a consequence of the wave packet ap- 
proach the particles are localized within a certain mo- 
mentum space volume. Consequently, the wavefunction 
decays into clusters which do not substantially overlap, 
and therefore the symmetrization has to be carried out 
only clusterwise. This is described in detail in Sect. 3.2. 
In Sect. 3.3 we demonstrate that clustering alone is 
not sufficient to reduce the numerical complexity in case 
of high momentum space densities. By means of simula- 
tions we investigate the impact of further factorizations 
on the correlation data. 



3.1 Metropolis procedure 

The method to generate an ensemble of test-particles 
with the statistical properties described in Sect. 2.2 is 
the standard Metropolis procedure [20]. At the begin- 
ning of each event only the phase space distribution Eq. 
(7) of wavepacket centres is defined. We now build a set 
of n test-particles in momentum space. Although their 
initial momenta are arbitrary, some specific choices may 
accelerate the speed of convergence to the correct statis- 
tics. We choose the initial momenta in a way which is 
dictated by the single-particle wave functions, i. e. the 
coordinates are sampled within the Gaussian probabil- 
ity distributions defined by the packets (j)i{p) Eq. (15). 
The advantage is that the ensemble of particles imme- 
diately satisfies the single-particle momentum spectrum 
Eq. (20). There are still, of course, no correlations cre- 
ated. These are introduced by the following procedure: 
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— Evaluate the value of the ra-particle density function 

l*(Pl,---,Pn)l^ 

— Do . . . 

— Change coordinate of one test-particle Pi — > Pi + 
^P 

— Evaluate density |^|new with the new set of coor- 
dinates 

— Accept new configuration with probability V = 

— . . . Loop over all particles 

After each particle was treated once, one sweep is fin- 
ished. It needs several sweeps until the coordinates of 
test-particles have converged to values which are com- 
patible with the momentum space probability distribu- 
tion |#|^. How many sweeps are needed depends on the 
size n of the system as well as on the specific properties 
of the density function. We found by "trial and error" 
that the rule 



N. 



sweep 



30 X 7x 



(37) 



yields results that do not change by further increasing 
the value of .ATs^eep- A second parameter of the simula- 
tion is the size of displacement ^p. It is a well known 
rule [21] that these displacements must not exceed the 
natural fluctuations inherent in the system. Otherwise 
unphysical states have a large probability to be reached 
and the acceptance rate of the Metropolis selection be- 
comes too small. The fluctuations are given by the widths 
(7p of the wavepackets in momentum space, and therefore 
we sample ^p within a uniform sphere of radius (7p. 

A series of simulations with different wavepacket 
widths (7o was carried out. Each simulation contained 
150000 events of 5 pions. The phase space distribution of 
the wavepackets was chosen by the simple parametriza- 
tions Eq. (4) and Eq. (19), with temperature T = 50 
MeV and source radius R = S fm. The results are dis- 
played in Fig. 2. 

The effective temperatures are in qualitative agree- 
ment with the calculation of the single-particle spectrum 
Eq. (20), but systematically smaller. Such a reduction 
is not observed in the single-pion system. It therefore 
seems that the interference phenomena present in the 5- 
pion system tend to decrease the pion energies relative to 
the expected single-particle spectrum. The effect is small 
(sa 4% for (Tq = 1.0 fm) but significant and becomes more 
pronounced the smaller the wavepacket width is in con- 
figuration space, i. e. the larger the overlap in momentum 
space. 

The correlation function is obtained by an eventwise 
evaluation of all pairs pj — pj . They contribute to the 
nominator of the correlation function. The denominator, 
that represents the uncorrelated 2-particle distribution, 
is obtained using event-mixing techniques. The standard 
correlation function Eq. (6) is fitted to the resulting dis- 
tribution, yielding an effective radius iZeff and a chaotic- 
ity A. For comparison, the correlation function Eq. (2) 
was solved numerically using the states Eq. (24) and Eq. 
(20) for a large number of wavepacket widths. Again, the 
standard correlation function was used to obtain radius 
and chaoticity parameter for each solution, yielding the 
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Fig. 2. The dependence of effective temperature Tgff (upper 
panel), effective radius i?eff (niiddle panel) and chaoticity A (lower 
panel) on the wavepacket widths CTq . The circles represent the sim- 
ulation results. Upper panel: The solid curve represents Eq. (21), 
the dashed line denotes the input temperature. Middle panel: The 
solid curve is in accordance to the exact solution of Eq. (24), the 
dotted curve is the approximation Eq. (29), the dashed line denotes 
the source radius. Lower panel: The solid curve is in accordance to 
the exact solution of Eq. (24), the dotted curve is the approxima- 
tion Eq. (33) and the dashed line denotes the constant chaoticity 
obtained in the approximation Eq. (28). 
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Fig. 3. The correlation function obtained by numerical (Monte 
Carlo) integration for the system used in Sect. 3 and <To = 1.8 
fm (circles). The solid curve is the fit of the Gaussian standard 
correlation function Eq. (6). 
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solid curve in the middle panel of Fig. 2. It is worth not- 
ing that the exact correlation function is not Gaussian, 
although the source density was chosen to be Gaussian. 
Figure 3 displays the correlation function corresponding 
to (7o = 1.8 fm (circles) and the fit of the standard cor- 
relation function (solid curve). Clearly, the Gaussian fit- 
function is not the best choice, but in order to keep the 
results comparable with the experimental data we de- 
cided to avoid any modifications of the common analysis 
procedures. The difficulties in the fit procedures gener- 
ate systematical errors, that in our simulations are larger 
than the statistical errors. In the figures, however, only 
the statistical errors are displayed. 

The obtained radii are larger than the one expected 
using the zero'th order approximation (Eq. (29), dotted 
curve in the middle panel of Fig. 2), but in agreement 
with the exact solution of the correlation function (solid 
curve). It is interesting that the exact results of the ef- 
fective radii can be parametrized with a function of the 
same type as Eq. (29), namely 



R. 



eff ■ 



IR^ 



Teff 



(38) 



yielding k= 3.05 ±0.01. 

In accordance to the discussions presented in Sect. 
2.2.2, the chaoticity A is not constant but decreases with 
increasing wavepacket width (lower panel). The solid 
curve represents the exact numerical solution, the dotted 
curve is the approximation Eq. (33) obtained by decou- 
pling momentum- and configuration space. The accuracy 
of the approximation is expected to decrease in case of a 
strong coupling of momentum- and configuration space, 
i. e. flow. Quantitative estimates have to be obtained by 
extending the simulations. 

The simulations were performed for a small system of 
only 5 pions. The generation of 150000 events required 
about 3 hours of CPU-time on a DEC 250^/^^^ worksta- 
tion. For larger systems the CPU-time increases dramat- 
ically. One way to overcome these problems is to split the 
system up into smaller subsystems. This is described in 
the following section. 



3.2 Cluster identification 

The direct evaluation of the ra-particle permanent Eq. 
(18) requires to sum up n\ terms. There are, of course, 
more intelligent algorithms available [22], but they still 
exhibit an exponential dependence of the complexity on 
the system size. Whereas in the case of determinants so 
called elimination algorithms like the Gaussian elimina- 
tion method with a complexity ~ n^ are available, a cor- 
responding procedure for permanents seems to be out of 
reach [23]. Faced with the problem to evaluate systems 
of much more than 10 pions, we have to find a way to 
circumvent the necessity of computing the complete n- 
particle permanent. The solution can be found in terms 
of an effective factorization. 

In order to understand this phenomenon, we consider 
the example of two particles. The probability density in 



momentum space is given by Eq. (25). Obviously the 
state factorizes if the second term |^2(pi)|^ |^i(P2)|^ and 
the interference term I vanish. Inserting Eq. (15) in Eq. 
(26) yields 



Re\I} = exp 



i^""' 



Pl)' + (P2-P2)' 



(P2-Pl)' + (P1-P2)') 



xcos((Xi-X2)-(p2-pi)/ft) . (39) 

If the distance of the wavepackets |P2 — Pi| is much 
larger than the width (7p, then necessarily one of the 
cross terms (Pj — Pj)^ aquires a value much larger than 
l/4(7p so that the exponential term and consequently the 
interference vanishes. The same argument holds for the 
term |^2(pi)|^ |^i(P2)|^- Thus the states are separated 
and a symmetrization does not influence the probability 
density Eq. (25). 

We want to exploit this result and split a given 
momentum space wavefunction into smaller parts. This 
means that clusters of overlapping wavepackets have to 
be identified. We can take advantage of the fact that 
the wavefunction in momentum space is stationary, c. f. 
Eq. (15), which means that this procedure has to be ap- 
plied only once per event. First of all one has to define 
under which conditions two wavepackets factorize. We 
have found that reasonable results are obtained using a 
cut-off value of 



|P2-Pl|<2.5<7p 



(40) 



which implies that wavepackets with distances smaller 
than 2.5 (7p are defined to overlap, otherwise they factor- 
ize. The second step is to construct a matrix [tableau) 
which indicates whether two states overlap or not. Fi- 
nally a tracer is launched to find a way through the 
tableau in order to locate the overlapping states. 

If the momentum space density of the system in- 
creases, it may happen that the wavepackets form a very 
large cluster [percolation). Taking into account the fact 
that the symmetrization of clusters with sizes larger than 
sa 10 takes too much time, the cluster has to be split up 
into subclusters by brute force, i. e. the wavefunction is 
forced to factorize. The influence of neglecting some of 
the interference terms on the correlation data is investi- 
gated in the following section. 



3.3 Unauthorized factorization of wavefunctions 

In order to keep the CPU-time within a tolerable range, 
it is necessary to factorize the wavefunctions of systems 
with very high momentum space densities, without tak- 
ing care whether the criterion Eq. (40) is satisfied or not. 
To examine the influence of such an operation on the 
correlation function, a system of 30 pions was simulated 
under different conditions: The maximum size Nmaic of 
the clusters was varied between 2 and 10. This means: 
If the cluster detection algorithm has identified a clus- 
ter of multiplicity n, then, if tx > Nmax., the cluster is 
divided into smaller subclusters. The effect of such an 
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Fig. 4. The distribution of cluster multiplicities identified by the 
cluster-tracer. Only a subset of them (here: clusters of size < 10, 
shadowed) can be treated in full symmetry, the rest has to be split 
up into smaller parts. 
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Fig. 5. Simulation results (triangles) for different maximum clus- 
ter sizes iVmax* Upper panel: The effective correlation radii. The 
solid line denotes the result of the numerical evaluation of Eq. (24), 
the dashed line is the source radius. Middle panel: The chaoticity 
A. Again, the solid line denotes the result of Eq. (24). Lower panel: 
The CPU-time per sweep. 



unauthorized factorization is of interest since most of the 
available codes generate pion correlations only in terms 
of 2-particle correlations Eq. (3). It is obvious that these 
methods become inaccurate if a certain phase space den- 
sity is exceeded and correlations of higher order than 2 
begin to dominate. 

In our simulation, the source was parametrized using 
the standard Gaussian Eq. (4) with radius R= h fm, the 
momentum distribution was the relativistic Maxwellian 



Fig. 6. Simulationresults (triangles) for different system multiplic- 
ities. Upper panel: The correlation radii. The dashed line denotes 
the input radius. Lower panel: The chaoticity A. The dashed line 
denotes the expected value A = 1, corresponding to the chaotic 



V{V) ~ exp 



^/V^ 



(41) 



with temperature T = 50 MeV. The wavepacket widths 
were (Tq = 1.8 fm. Figure 4 shows the distribution of clus- 
ter multiplicities. A large amount of clusters with sizes of 
around 25 particles was found, indicating that percola- 
tion phenomena play a significant role. Figure 5 displays 
the effects appearing if the clusters are broken into pieces 
of .ATniax = 2, 3, ...,10. The upper panel shows the effec- 
tive radii obtained in the same way as described in Sect. 
3.1. The forced factorization obviously has no influence 
on the extracted radii. These also agree with the numer- 
ical solution of the 2-particle correlation function Eq. 
(24), which yields R^s = 5.70 ± 0.05 fm (solid line). In 
the middle panel a strong dependence of the chaoticity 
parameter A on Nmaic is displayed (triangles). The result 
of the numerical solution of Eq. (24) (A = 0.94 ± 0.01, 
solid line) is not reached by any simulation. This implies 
that the neglection of interference terms in the simula- 
tion procedure makes an accurate determination of the 
chaoticity (and therefore of the number of degrees of free- 
dom in configuration space) impossible. This has to be 
taken into account when only 2-particle correlations are 
used to simulate events. 

There exists an additional effect which is described 
by Zajc [9]: In case of high phase space densities also 
the analysis suffers from the fact that the 2-particle cor- 
relation function ignores the contributions of higher or- 
der correlations. This phenomenon leads to a reduced 
value for A and the effective radius. In order to study 
the influence quantitatively, a series of simulations with 
increasing system size was carried out using the quan- 
tum statistical code presented in [8]. This code is based 
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on the plane wave approximation (35), i. e. it per defini- 
tion generates pions from a chaotic source and yields the 
correlation radius iZeff = R- On the other hand, it in- 
cludes Bose-symmetrization up to all orders. Thus, this 
method offers the best choice to study the effects which 
appear exclusively in the analysis procedures. In the sim- 
ulations, the same phase space parametrization as above 
[R = 5 fm, T = 50 MeV) was used, the results are 
presented in Fig. 6. The fact that the radius comes out 
somewhat too large for all simulations is a consequence 
of the discretization of the momentum space into cells. In 
agreement with the predictions of Zajc, the radius as well 
as A decrease with increasing system size, for the latter 
the effect is more pronounced. For a system of 30 pions 
this is already visible and therefore the extracted A val- 
ues in Fig. 5 are expected to saturate at approximately 
A(^max = 30) sa 0.75 instead of 0.94. We want to em- 
phasize that here, in contrast to the simulations shown in 
Fig. 6, the deviation of A from 1 is due to a combination 
of three effects, (1) the omission of interference terms in 
the simulation, (2) the omission of higher order corre- 
lations in the analysis and (3) partial coherence, where 
the contributions from (2) and (3) are independent on 
-^max- The contribution from (2) also depends on the 
phase space volume, because the number of momentum 
space cells is proportional to the product R'^T [8]. The 
larger the number of cells, the smaller the probability 
to find more than two pions in identical cells, which im- 
plies that the amount of reduction becomes smaller. For 
a system similar to the 200 AGeV Pb + Pb experiment in 
CERN (e. g. iZ = 6 fm, T = 150 MeV, pion multiplicity: 
1000), we obtain R = (5.6 ± 0.1) fm and A = 0.75 ± 0.04 
for a chaotic source. This implies that for an accurate 
measurement of the chaoticity higher order correlation 
functions have to be included. 

Finally, the lower panel in Fig. 5 displays the needed 
CPU-time per sweep. It is clear that the evaluation of 
clusters with Nmaic much larger than 10 is prohibited by 
the exponential growth of the operations to be carried 
out by the algorithm. The generation of 2000 events un- 
der NmaiL = 10 - conditions needed about 15 hours on a 
DEC 250^/266 workstation. 

In this and the foregoing section two ways of reducing 
the calculational effort in the simulation were presented: 

1. Truncation of the ra-particle wavefunction into clus- 
ters in accordance with the criterion of effective fac- 
torization. This operation does not effect the corre- 
lation properties. 

2. A split up of the clusters into subclusters of sizes 
< -^max- This operation is unauthorized in the sense 
that interference terms are lost. Although the cor- 
rect source radii can be extracted even in case of 
-^max = 2, the chaoticity A in the correlation function 
decreases dramatically. 

In addition, the factorization of the system by analyzing 
only 2-particle correlations also effects A in a way that 
an accurate measurement of the chaoticity under e. g. 
CERN conditions becomes possible only when higher or- 
der correlation functions are used. 



4 Molecular dynamic method 

The Monte Carlo procedure presented above allows to 
study a large class of physical scenarios, including col- 
lective flow and dynamic changes of the source. However, 
these situations are correctly treated only when the pi- 
ons can be regarded as free, which implies that there are 
no forces acting on them after freeze-out. This is due to 
the fact that the Monte Carlo procedure is not capable 
of providing a physical timescale, instead it is designed 
to generate statistical ensembles under the condition of 
stationarity. 

On the other hand, final state interactions of long 
range can occur in heavy ion collisions, for example 

1. the Coulomb interaction between pions and the cen- 
tral fireball, 

2. the Coulomb interaction between pion and pion, 

3. the Coulomb interaction between pions and the spec- 
tators. 

These interactions require a time dependent treatment 
in configuration space, i. e. an integration of the equa- 
tion of motion. This is the domain of molecular dynamic 
procedures. 

At this point we are confronted with a fundamental 
problem: Quantum "mechanics" does not provide us with 
an equation of motion. Instead it contains a number of 
rules for the time development of probability densities. 
Thus the evolution of a system has to be found by a 
numerical solution of the ra-particle Schrodinger equa- 
tion, which is very time consuming and inadequate to 
our goal. It is remarkable that there exists an extension 
of the usual quantum mechanical formalism [24] (for a 
thorough treatment including applications see [26]) that 
allows us to formulate rules for the motion of certain test- 
particles quite similar to the particles used in the Monte 
Carlo method. In the following section the application 
of this formalism to our free-particle model (Sect. 2) is 
demonstrated. Sect. 4.2 contains an approximative treat- 
ment of Coulomb interaction between pions and fireball. 



4-1 B ohm's quantum theory of motion 

The system in configuration space is described by the 
7x-particle wave function Eq. (13). We consider a set of 
n test-particles with phase space coordinates 

{(Xi(t), pi(t)), (X2(t), P2(t)), . . . , (X„(t), pnit))} . (42) 

Further we may choose a set of initial conditions for the 
configuration space coordinates of the test-particles 

7'(X1,X2,...,X„,0)= |!?(xi,X2,...,x„,0)|2 (43) 

as well as for the momentum space coordinates 

Pj(0) = Vj5(xi,...,x„,0), (44) 

where 



2i ^\W* 



(45) 
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is the phase of the wavefunction and j = 1,2, ... ,n. If 
the wavepackets are small enough so that their overlap 
can be neglected, the initial distribution reduces to 



7'(x,-,0)=|V-,-(x,-,0)| = 

and 

Pi(0) = Pi • 



(46) 



(47) 

It was shown by Bohm [25] that application of the fol- 
lowing equation of motion to the test particles yields in 
the limit t — > co a phase space distribution which coin- 
cides with the configuration space and momentum space 
density Eq. (18): 



d'^jjt) 



dt^ 



-VjQ(xi,X2,...,X„,t) 



with 



k = l 



2m 



R 



(48) 



(49) 



R = (^*^y/'^ denotes the amplitude of the wavefunc- 
tion. Obviously, Eq. (48) is an equation of motion of 
Newton type, but the potential Q, called quantum po- 
tential, is nonlocal since it contains the (symmetrized) 
amplitude. We want to emphasize that 

1. the equation of motion Eq. (48) is not based on a 
semiclassical approximation but is of pure quantum 
mechanical nature, 

2. Bohm's picture of quantum mechanics leads to the 
same observable predictions for the test-particles as 
the conventional approach for the results of measure- 
ments. Thus all ensemble properties obtained in Sect. 
2.2 remain valid for the test-particles. 

The quantum potential Eq. (49) for the wavepacket 
model can be expressed analytically but exhibits a quite 
complicated structure for many-particle systems. Here 
we only want to discuss the most simple case of one 
particle. Considering a wavepacket Eq. (11) at rest (P = 
0) and centred at X = we obtain 



Q(x,t) 



h^ 



irruT^ 



'it) 



20-2 (t) 



(50) 



where a[t) = ■^s[t)s*[t) is the time dependent ampli- 
tude of the complex width s[t) Eq. (12). If we are inter- 
ested in the average initial quantum potential of a large 
ensemble of single-particle systems, we have to integrate 
over the wavepacket-density distribution 



(27r<7o) 



-3/2 



exp 



20-2 



|V-(x,t = 0)|2 

obtaining 

(Q(t = 0)) = y"Q(x,0)|V-(x,0)|2d^ 

and comparing with Eq. (21) one finds 



3fe2 
8m(72 



(51) 



(52) 



(53) 
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Fig. 7. Contour plots of the quantum potential corresponding to 
a 3-particle system. Beginning with the upper right panel and pro- 
ceeding clockwise, the pictures display the potential after 0, 10, 15 
and 25 fm/c. Two wavepackets were fixed in the ajy-plane, a third 
was used to scan the plane. The labels are in fm-units. 



This implies that the quantum potential in Bohm's pic- 
ture of individual particles manifests itself (in the en- 
semble picture) in terms of the quantum temperature 
we have found to be a consequence of Heisenberg's un- 
certainty in momentum space. How this determines the 
dynamics of the particles can be seen by looking at the 
quantum force 



VQ 



ft2 



47710-4 (t) 



X(t). 



(54) 



If the particle is placed at the centre of the wavepacket 
(x(0) = 0), the quantum force disappears and the par- 
ticle follows a classical trajectory (Ehrenfest's theorem). 
Otherwise, the force always points away from the wave 
packet centre, i. e. in a large ensemble the test-particles 
are accelerated isotropically in all directions, indicating 
the dispersion of the wave packet. In addition, the force 
decreases rapidly with increasing spreading (o""'* depen- 
dence). What we can learn from this example is that 
the quantum force accounts for Heisenberg's uncertainty 
principle, but from an individual rather than from an 
ensemble average point of view. 

To get an impression. Fig. 7 shows the quantum 
potential of a 3-particle system and it's time develop- 
ment. To obtain the pictures, two wavepackets with ini- 
tial widths (7o = 2.0 fm and centred test-particles were 
fixed in the asy-plane. Since they are at rest, the only dy- 
namics visible is due to dispersion. A third wavepacket 
(with centred test-particle) is used to scan the asy-plane 
at different times. Therefore, the pictures display the po- 
tential seen by the test-particle of the third wavepacket. 

For numerical simulations, there exists a different and 
less expensive way to solve the equation of motion: The 
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velocities of the particles can be obtained by Eq. (44), 
which we can rewrite as 



Vi(<) 



h 



m\W\ 



Im{^*Vj^} . 



(55) 



The force acting on the particle in the discrete time inter- 
val dt = tn — tn-i can be expressed as — VQ = mdw /dt 
with dv = v(t„) — v(t„_i). The gradient of the perma- 
nent !? is obtained in the following way: 

Vfc!? = Vj, I ^V'<r(l)(xi)V'<^(2)(x2)---V'<^(„)(x„) J 

= X] ^*(V'<r(fc)(xfc)) V'<^(l)(xi) • • • V'<^(fc-l)(xj,_i) 
xV'<7(fc + l)(xfc + l) • • •V'<r(n)(Xn) 

i = l \ a^i 



n 

^V,(V-,(x,))!?('''=)(xi,...,x„) 



(56) 



where !?(''*''(xi, . . . , x„) is the permanent of the ma- 
trix after canceling the i'th row and fe'th column in 
!?(xi, . . . , x„). We note that gradients of permanents are 
sums of permanents of the submatrices. For the compu- 
tation of the permanents the code presented in [22] is 
used. 

It is clear that the presence of such irregular poten- 
tials as displayed in Fig. 7 requires very small timesteps 
for the integration of the equation of motion. This sim- 
ulation method therefore proved to be usable only for 
small systems of less than sa 10 particles. However, 
there is some kind of separability inherent in the sys- 
tem: Figure 7 suggests that the quantum potentials of 
the single wavepackets disturb each other only in case 
of overlap, and since the packets move away from each 
other in a realistic simulation, it is possible to divide 
the wavefunction into subparts in the same way as de- 
scribed in Sect. 3.2. The separation procedure has to be 
repeated frequently, since the configuration changes with 
the progress of time. 



^.2 Simulation including Coulomb interaction 

It was our goal to include final state interactions into 
the simulation, and we present the example of the pion- 
fireball Coulomb interaction. For this purpose we de- 
scribe the charge distribution by a homogeneous sphere 
of rms-radius r = 8.6 fm and total charge Z = 100, cor- 
responding to 100 positive nucleons in the participant 
region. For the equation of motion, we use the following 
approximation : 






-V,(Q + T^) 



(57) 



where V denotes the classical (static) Coulomb potential. 
The wavepackets are also exposed to the Coulomb field, 
because they represent a cloud of charge Z = 1, but 
their shapes remain unefFected. This picture is correct 
only as far as the Coulomb field is homogeneous on the 
scale of the wavepacket diameter. Otherwise the quan- 
tum potential becomes a function of V. In this sense, Q 
in Eq. (57) can be regarded as the zero'th order Taylor 
approximation of Q{V). 

For the simulations, wavepackets with initial widths 
(To = 1.8 fm were used. Again, the phase space distrib- 
ution of the wavepackets was obtained by the parame- 
trizations Eq. (4) with radius R = 7 fm (in accordance to 
the rms-radius ?■ = 8.6 fm of the charged sphere) and Eq. 
(41) with temperature T = 50 MeV. The initial phase 
space distributions of the test-particles were chosen in 
accordance to conditions Eq. (43) and Eq. (47), which is 
simple as long as the configuration space density is suffi- 
ciently small that the overlap of the initial wavepackets 
can be neglected. If this is not the case, a Monte Carlo 
procedure similar to the one discussed in Sect. 3 has to 
be applied to obtain the initial configuration space dis- 
tribution Eq. (43). 

In the simulations, 10^ events for neutral, positive 
and negative pions were generated. The adaptive time 
step was chosen in a way that the energy gain (resp. 
loss) of one particle did not exceed 1 MeV per timestep. 
Pion-pion interactions were not taken into account, but 
can be, of course, included in a similar way as the pion- 
fireball interaction. Each event contained a number of 6 
pions, the CPU-time needed was about 1 s per event on 
a DEC 250^/266 workstation. 

While the quantum forces can aquire quite high val- 
ues (up to several hundreds of MeV/fm) during the first 
phase of the evolution, they tend to decline rapidly in 
the course of the wavepacket dispersion (see Eq. (54) and 
Fig. 7), so that there exists a moment when they vanish 
and the equation of motion becomes classical. We found 
that a cut-off value of |F| = 0.01 MeV is a sensible choice 
to stop the further evaluation of the event, since from 
that point on the momentum space distribution of the 
test-particles does not change substantially. For the neu- 
tral particles, an average number of 580 timesteps was 
needed to reach this point (which we may call "quan- 
tum freeze-out"), in terms of the system time this value 
corresponds to sa 130 fm/c, but both values can fluctu- 
ate substantially from event to event. In case of charged 
pions, these values are higher due to the long range char- 
acter of the Coulomb force. 

Figure 8 displays the resulting normalized momen- 
tum spectra. It is not surprising that the spectra of the 
TT"'" [w~ ) are somewhat enhanced (suppressed) for low 
momenta in comparison to neutral pions, a consequence 
of the energy gain (loss) due to the central Coulomb 
potential. The results of the correlation analysis are dis- 
played in Table 1. The values were obtained by fitting the 
standard correlation function Eq. (6) to the correlation 
signal Eq. (2) extracted from the simulated data. Within 
the errors, the w° - radii are in agreement with the re- 
sults obtained by the numerical evaluation of Eq. (24), 
labeled as 7r°(th. ) in Table 1. The w~ data, however. 
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Table 1. Extracted radii i?efF and chaoticity parameter A for neu- 
tral pions (7r°), positive charged (tt""") and negative charged pions 
(tt") and for different momentum cuts. For comparison, results 
of the numerical evaluation of Eq. (24) are displayed in the last 
column. 



p (MeV) 




7r° 


■K+ 


TT 


7r° (th. ) 


0- 600 


^eff 

A 


7.7±0.3 
l.OiO.l 


7.3±0.3 
l.OiO.l 


9.1±0.3 
0.9±0.1 


7.44±0.04 
0.96±0.01 


0- 120 


^eff 
A 


7.3±0.5 
l.OiO.l 


8.1±0.7 
1.1±0.2 


10.8±0.8 
1.0±0.2 


7.35±0.04 
0.96±0.01 


120- 200 


^eff 
A 


7.5±0.5 
0.9±0.1 


7.1±0.4 
0.9±0.1 


8.2±0.6 
0.8±0.1 


7.43±0.04 
0.96±0.01 


200 - 600 


^eff 
A 


7.0±0.6 
1.0±0.2 


6.4±0.5 
0.8±0.1 


7.5±1.0 
0.9±0.2 


7.44±0.04 
0.96±0.01 
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Fig. 8. Normalized momentum spectra for 7r° (circles), ir'^ 
(squares) and 7r~ (triangles). 



exhibit a significant radius increase for low momentum 
pions and approach the value of the nominal radius for 
high momenta. Integrated over all momentum slices, the 
w~ - radius is about 20 % larger than the w° - radius. 
On the other hand, the tt"'" - radii appear slightly en- 
hanced in the low momentum range and drop below the 
w° values in the high momentum range. The integrated 
radius agrees, within the errors, with the integrated w° - 
radius. Within the statistical errors, the Coulomb inter- 
action between charged pions and the fireball does not 
appear to effect the chaoticity parameter A significantly. 
It is an open question whether or not a more accurate 
treatment of the quantum potential, i. e. the inclusion of 
higher order Taylor terms in Eq. (57), may change the 
results presented in Table 1. But they may indicate the 
direction in which the deduced radii will change when the 
Coulomb interaction between the pions and the source 
is no longer negligible. They are compatible with the 
results obtained by Barz [27], who has solved the Klein- 
Gordon equation with Coulomb potential numerically by 
using a partial wave expansion technique. Bohm's the- 



ory provides an alternative way to approach the quantum 
mechanics of a system of particles with residual interac- 
tions: Approximative methods like perturbation theory 
can be applied to the quantum potential, which deter- 
mines the dynamic changes of the system on a funda- 
mental level. 



5 Conclusion 

In this work we have developed a method which trans- 
forms a given classical distribution of particles in phase 
space into a corresponding Bose-Einstein quantal distri- 
bution. In Sect. 2 we have demonstrated how this can 
be achieved by first building a symmetrical state from 
the original particle distribution and then introducing 
a set of test particles as a representative of the quan- 
tum mechanical ensemble. We have shown in Sect. 2.2 
how the statistical properties of the resulting mixed en- 
semble deviate from the original classical ensemble: As a 
consequence of Heisenberg's uncertainty principle single 
particle distributions in configuration- and momentum 
space are smeared out. The relevant parameter is the 
initial wavepacket width (Tq, which represents the degree 
of information loss one is willing to accept either in con- 
figuration space or momentum space. 

The two-particle correlation function for wavepackets 
is complicated and has to be treated numerically, but we 
have demonstrated that even for the simple example of 
a static, Gaussian and thermal source 

1. the correlation function is not Gaussian (Eq. (24) and 
Fig. 3), 

2. the correlation radius disagrees with the geometrical 
radius (Eq. (29)), 

3. and the chaoticity A is a complicated function of the 
system parameters, in particular it depends on the 
extensions of the wavepackets relative to the source 
size (Eq. (33) and Fig. 2). 

In Sect. 2.3 we have shown that the usual plane wave ap- 
proach discouples the configuration space from the mo- 
mentum space. It violates Heisenberg's uncertainty prin- 
ciple and does not provide a physical explanation for the 
chaoticity A. In contrast, the wavepacket picture allows 
to relate the observable chaoticity with the localization 
of the pions in the source, i. e. their mean free path. This 
implies that A may contain information about the cross 
section for pion absorption in nuclear matter. 

In Sect. 3 we presented a Monte Carlo procedure to 
treat systems of free (noninteracting) particles in mo- 
mentum space. The consistency of the results with the 
expected ensemble properties was tested in Sect. 3.1. One 
interesting result is that the temperature deduced from 
the measured momentum spectrum of pions is smaller 
than the theoretical temperature of the single particle 
spectrum. A plausible reason would be the reduction of 
the degrees of freedom caused by the interference. 

The large numerical effort needed to evaluate sys- 
tems of large pion numbers motivated us to develop cri- 
teria under which the system may separate. In case of ef- 
fective factorization, the correlation phenomena remain 
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undisturbed (Sect. 3.2). However, the study of high mo- 
mentum space densities requests a further, unauthorized 
factorization of the states. We have shown that also this 
factorization does not lead to a noticeable change of the 
effective correlation radii, but the chaoticity A becomes 
smaller as a consequence of omitting interference terms 
(Sect. 3.3). Therefore we argued that usual codes which 
only use pion pairs to generate correlated states are not 
capable of treating the chaoticity in a consistent way. 
We further want to point out that our approach to para- 
metrize the system in terms of the standard pion source 
and the Maxwell distribution and to use the standard 
correlation function is an oversimplification, since it con- 
tains no more physics than "radius" , "temperature" and 
"chaoticity". Phenomena like flow, time dependencies of 
the source, resonance decay and final state interactions 
should be included since they are readily observed in 
heavy ion collisions. The model we have presented is 
capable to contain these cases since it can be applied 
to any phase space distribution. An adequate analysis 
of these phenomena, however, may ask for finer observ- 
ables than just 2-particle correlations. As an example we 
have demonstrated that an accurate measurement of the 
chaoticity A under conditions realized in the CERN Pb 
+ Pb (200 AGeV) experiment, i. e. multiplicities in the 
order of 1000 pions, requires the analysis of higher order 
correlation functions. 

Interactions of long range such as the Coulomb in- 
teraction can not be treated exclusively in momentum 
space. Instead, a time dependent calculation in config- 
uration space in the framework of molecular dynamics 
has to be performed. The extension of the usual quan- 
tum mechanics given by Bohm offers a method to pro- 
ceed from the laws of ensemble averages to the laws of 
a single representative (Sect. 4.1). The quantum forces, 
however, are of quite irregular nature and therefore a nu- 
merical treatment will usually be restricted to systems 
of low multiplicities. 

The treatment of pion-fireball Coulomb interaction 
in Sect. 4.2 should be regarded as a first attempt to han- 
dle final state interactions without neglecting the quan- 
tum mechanical features of the system. The resulting 
discrepancy between the correlation radii of tt"'" respec- 
tively 7r~ pairs, which is largest for low pion momenta, 
has been observed in experiments [19]. Our simulations 
indicate that the effects on the radii are not of opposite 
size for the two pion charges, i. e. taking the average of 
both values as estimate for the source radius is not an 
adequate choice. 



the statistical weight of the i'th cell. We want to show 
that Vc = "^liK + 1), independent on N. 

Be N = 1. The occupied j'th cell has the statistical 
weight w(j) = 2/[K-\-l), whereas all other cells have the 
weights w(i) = 1/{K + 1), i ^^ j. For the second pion, 
we therefore obtain 



V,{N = 2) = w{j) 



K + 1 



(59) 



We now assume that Eq. (59) holds for a certain N > 2, 
which implies that 

, , m^ yj, n(i)(n(i) - l)/2 2 , , 

VAN) = —^ = ^'-^ } '^ ^ { , — '-^ = . (60) 

^ ^ m N{N -l)/2 K + l ^ ' 

We now consider the N + I'th pion and calculate how 
many coherent pairs it is expected to build with the other 
pions: 



K 



K 



(me)new = ^ ^{i) n{i) = ^ — 



{n{i) + l)n{i) 



i = l 

We therefore obtain 



i = l 



N 



(61) 



V^^N + 1) = "^- + <"^-^- 

^K ( n(i)(n(i)-l) (n(i) + l)n(i) \ 
Z^i = l \ 2 ~'~ K+N ) 



N{N-1) (JV+1) 
2 (JV-1) 



l^i = l 2 y-^ K+N ) ^ K+N 



N{N-1) (JV+1) 
2 (JV-1) 

K 



,(62) 



where we have used that X^^-^ ra(i) = N . If we now sub- 
stitute Eq. (60), we get 

VAN^l) = ^^-U' ' ^(^-^) 



(63) 



K + l\ K + Nj {N + 1) 
4 



{K + N){N + 1) 
2 



K + 1 
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6 Appendix 

Be K the number of cells, N the number of pions which 
are already distributed into the cells, m the number of 
pairs, and rric the number of coherent pairs, i. e. the num- 
ber of pairs that can be built from pions which occupy 
identical cells. We define Vc as probability to find a co- 
herent pair, n[i) the number of pions in the i'th cell and 

n{i) + 1 



,{i) 



K + N 



(58) 
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